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Abstract 

We consider three-dimensional clusters of equal-volume bubbles packed around a central 
bubble and calculate their energy and optimal shape. We obtain the surface area and bubble 
pressures to improve on existing growth laws for three-dimensional bubble clusters. We dis- 
cuss the possible number of bubbles that can be packed around a central one: the "kissing 
problem", here adapted to deformable objects. 

Pacs numbers: 82.70.Rr, 83.80.Iz 

1 Introduction 

1.1 Motivation 

Bubbles, such as soap bubbles, are objects with simple geometry and physical properties. But 
when two or more bubbles cluster together, how well do we really understand their properties? 

The limiting case of a cluster of many bubbles, known as a foam, is usually approached with 
continuum approximations. An understanding of foam properties such as aging, due to gas dif- 
fusion, and structure is a problem of fundamental interest stimulated by the need to predict the 
behaviour of foam in industrial applications. From carbonated drinks to the processes used to ex- 
tract gold ore from the earth, foams are an important part of our lives with various industrial uses 

fli. 

The alternative to the continuum description, described here, is an approach based upon the 
study of finite clusters of bubbles. Its advantage is the ease with which we can obtain precise 
structural information. A further benefit of studying finite, rather than infinite or periodic, foams 
is that the bubbles are not "frustrated", so that we get a measure of their free shape, rather than one 
influenced by long-distance correlations between bubbles. 

This has been demonstrated convincingly in two-dimensions (2D), where exact results exist for 
two problems of paramount interest: 
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• The Kelvin problem: what is the least energ y (equivalent to line-length) structure of equal- 
size bubbles that fills space? In 2D, Halej OD proved that this is the familiar honeycomb 
structure. In 3D, where the problem is one of minimizing surface energy or area, no such 
exact result exists. Kelvin [4] gave a candidate structure, still believed to be the best for a 
structure containing identical cells, although in the general case it has since been beaten by 
the Weaire-Phelan structure JH] consisting of bubbles of two different types. The important 
quantity in this problem is the surface area of each face of a bubble of unit volume, or 
equivalently the normalized total surface area 

• Growth laws: how does a foam age, or coarsen, due to gas diffusion across its surfaces? 
The 2D result, due to 0], says that the growth-rate of a bubble (of area A) is directly linked 
to its number of sides, n: dA/dt (n — 6). That is, it depends upon bubble topology only, 
irrespective of the precise geometry. In 3D, the growth law is written [Q]: 

3 dV 2 / 3 _ 1 ApjSj 

2D eff dt 2 i yV3 1 ' 

where the sum is taken over each face, which has a pressure difference Api and area Sj. F> e ff 
is an effective diffusion coefficient. Again, the normalized area appears to be important, but 
does the 3D growth law depend only on the bubble topology? In fact it does not, but it may 
make sense to express the average growth-rate of F-faced bubbles as a function of F only if 
the dispersion about such a law is small. 



1.2 State of the art 

The study of 3D foam coarsening was pioneered by Glazier [7], who used a 3D Potts model to 



numerically simulate foam coarsening. He proposed a linear growth law T^ApiSi oc (F — est) for 
bubbles with a number of faces F from 6 to 57 (and even from 4 to 60, with some numerical 
uncertainty). Similar linear laws were observed in subsequent experiments involving optical to- 
mography and reconstruction using the Surface Evolver [8] (F between 9 and 16), and magnetic 
resonance imaging experiments [9, 10] for F from 4 to 26. 

This growth law was refined by three detailed results presented by Hilgenfel dt et al. I Till : first, 
an approximate analytical formula based upon regular F-faced polyhedra with curved faces: 
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where %p = 2 tan 1 \/4sin 2 (7l/r|/7) — 1 and r\p = 6 — 12/F is the number of edges per face. For 



large F, this shows a square-root dependence, G(F ^> 1) = 2. 14yF — 7.79, effective for F greater 
than about 15. Second, a (non-explicit) correction for non-regular faces; third, nu merical (Surface 
Evolver) simulations for foams containing bubbles with F from 5 to 42. Recently, Cox and Fortes! 
JT2I1 also used the Surface Evolver to calculate numerically the structural properties of single "reg- 
ular" bubbles with surfaces of constant mean curvature; this gave information for certain values of 
F between 2 and 32. 



1.3 Outline of this paper 

Here, we study clusters consisting of one bubble surrounded by F others, each with prescribed 
volumes. This constitutes a finite cluster with free boundary conditions: this represents a realistic 
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foam surrounded by air, in contrast to the idealisation used to derive © Hilgenfeldt et al. 1 11]. 
We chose such an approach, which neglects long-distance correlations between bubbles, because 
it should provide more physical insight than existing experiments and simulations, and enable 
more precise calculations than the analytical approach. Within this "mean-field" choice, all results 
presented below are highly accurate, without approximations. Moreover, in principle we should 
have access to all physically realizable values of F. 

In the course of our study of the equal- volume case, we encountered what we call "the kissing 
problem for (deformable) bubbles". Our simulations allow us to ask: how many deformable (dry) 
bubbles can be packed around one other? The original kissing problem, discussed by Gregory and 
Newton, was: how many identical hard spheres can surround one other, each touching the central 
one LI 3] ? In two-dimensions the answer is obvious and well-known - only six hard discs can 
be packed around one other, in the familiar honeycomb arrangement. For the three-dimensional 
problem, consideration of the angle subtended by each sphere at the central one suggested that 
the maximum number could be as high as 14, but Newton was correct in believing that only 12 
neighbours are possible [14]. We will present arguments suggesting that for bubbles these critical 
numbers are 12 (2D) and 32 (3D). 

The plan of the paper is as follows. We first describe our method of cluster preparation and 
relaxation. There are limits, for each set of given bubble volumes, to the values of F for which 
stable clusters exists. In the equal-volume case we offer a solution to the kissing problem. We then 
analyse in more detail the shape and growth-rate of many different bubbles, and present predictions 
about coarsening and quantify the spread of the growth-rate about the growth law 0. 



2 Definitions and methods 

We take a central bubble of volume V c and surround it with F bubbles, each with the same volume 
V; this is the natural extension into 3D of the 2D "flower" of II 1511 . To create and equilibrate such 
a cluster, we use a Voronoi construction with VCS iflfill and then the Surface Evolver lfl7ll . as 
follows. 

We must first make a choice about the topology of the cluster. Since we wish to create the 
cluster using a Voronoi routine, we must first choose an arrangement of F + 1 points about which 
to create bubbles. 

We first place a point at the origin of a sphere of radius 1 . Then the Voronoi points are placed 



at the positions given by the solution of the "covering radius problem" [13]: the arrangement of 
F points on the unit sphere that minimizes the maximum distance of any point from its closest 
neighbour. Candidates to the solution of this problem have been given by [ 18] for F from 4 to 130, 
which is exactly what is required for our purpose. Note that this is not the only way to pack the 
F Voronoi points, but appears (partly with hindsight) to have been a good choice - it gives all the 
arrangements we know to expect, e.g. for F = 6, 12, 32. 

We truncate the Voronoi diagram by adding 3F points at a radial distance of 2 from the origin. 
We ensure that these outer points are at least a distance 2e/v / 3~F apart, decreasing £ from 1 until a 
solution is found, usually at around £ = 0.8. This data is put through the VCS software; the output 
file is then transferred to the Surface Evolver, version 2.18d. We use two levels of refinement and 
quadratic mode, to obtain a high level of accuracy - we estimate all values to be accurate to at least 
four decimal places. 

We compute the following quantities for the i th face (i = 1, • • • ,F) of the central bubble: its 
number of sides n ; , area Si and pressure difference Ap ; . Then for the whole bubble we record 
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Figure 1: Examples of the clusters considered here with V c = V = 1 in each case: (a) cluster of 
F = 13 outer bubbles (that is, a total of F + 1 = 14 bubbles) still attached; (b) the central bubble, 



drawn to a different scale, with F = 13 faces (this is a Matzke cell 12211 with one square, ten 
pentagonal and two hexagonal faces); (c) a bubble with F = 26 faces (all pentagonal or hexagonal) 
- note its departure from approximate sphericity, described in 33.21 

its volume V, its normalized total line length L/V 1 / 3 , its normalized surface area 

S/y2/3 (where 

S = T-iSi), and its growth-rate through eq. O, which we plot as a function of F. 



3 Topology and limits for equal-volume clusters 

We first consider the case where the volume of the central bubble is equal to that of its neighbours, 
V c = V. Examples of such monodisperse clusters are shown in figure[l]for F = 13 and 26. This 
illustrates that despite the rather symmetric initial condition (putting points on a sphere) we can 
still obtain significantly skewed bubbles after relaxation. 



3.1 The kissing problem for 2D bubbles 

For completeness, we consider first the two-dimensional problem. How many 2D bubbles can be 
packed around one other of the same area? 

Our initial pattern is that of the "flower" clusters introduced recently lfl5lH9[Eoll . It consists 
of a central cell of area A c surrounded by F identical petals of area A. A symmetric example with 
F = 12 petals and A c = 2A = 2 is shown in figure|2a). A priori, one could imagine that the number 
of petals could increase without limit, with the F sides of the central bubble becoming increasingly 
curved. 



However, Weaire et al. 1 15] showed that for F > 6 there is a "buckling" instability at a critical 



ratio of the bubble areas given approximately by 

A c ./A^0.04(F-6) 2 . 

For unit areas and F > 6 + (0.04) ~ 1//2 =11, the symmetric shape is therefore no longer stable, the 
flower becomes "floppy" and many modes of buckling, corresponding to different shaped central 
bubbles, are possible (all with the same energy). An example for F = 12, in which there is an 
elliptical mode of buckling, is shown in figure|2tb). 

Is it possible to pack even more bubbles? We find that for F > 12, any of the buckled configu- 
rations of clusters with unit areas are unstable to a topological change caused by the length of one 
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a) b) c) 




Figure 2: (a) A symmetrical 2D flower cluster with F = 12 petals and A c /A = 2. (b) One of the 
possible stable buckled states of the same cluster with A c = A. (c) One of the possible "ejected" 
states W20I1 with F = 13 petals and A c = A. 

of the internal edges shrinking to zero lEoll . An example is shown in figure Efc) for F = 13, for 
which three bubbles are "ejected" in an equilibrium configuration with A c /A = 1. 

Thus we conjecture that the maximum number of bubbles that can touch the central one is 12. 
This is twice the value for hard discs. 

3.2 The kissing problem for 3D bubbles 

In three dimensions the idea is the same. In principle one could imagine that there should be no 
limit to the number of bubbles which will fit around the central one, albeit with the latter being 
hugely distorted. However, since the area of each of the five-sided faces shrinks as F increases, 
our simulations of bubbles with unit volumes, V c /V = 1 do not find a stable cluster for all possible 
values of F. In fact, we could only find clusters for 5 < F < 32. That is, we cannot obtain a bubble 
with more than 32 faces and volume equal to that of its neighbours which satisfies Plateau's laws 
after energy (surface area) minimization. 

For most values of F, the shrinkage of five-sided faces is accelerated by an ellipsoidal distortion 
of the central bubble (see figure [He)), due to the asymmetric location of the pentagonal faces 
amongst the hexagonal ones. Might there be a discontinuous buckling transition for 3D clusters? 
As a result of further simulations, we believe not: this asymmetry, and the consequent elliptical 
deformation of the central bubble, means that the transition to the asymmetric pattern is continuous. 

It is interesting to note that the case F = 32 is special: it is probably the most symmetric cluster 
for F > 12 - it corresponds to the C^o fullerene. Hence, by analogy, one might expect that stable 
clusters with unit volumes exist for higher order carbon structures. We tried Cso (F = 42) and the 
elliptical C72 (F = 40) and didn't find them to be stable. We thus conjecture that no more than 32 
bubbles can touch the central one: 32 appears to be the "kissing" number for 3D bubbles. Recall 
that for hard spheres the kissing number is 12. 

4 Shape, pressure and growth-rate 

We next analyse in detail the statistics of the bubbles found in our simulations. 
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Figure 3: The normalized surface area s/v 2 ' 3 varies, albeit over a small range, non-monotonically 
as F increases. Inset: zoom over the range F = 11 to 16. Data is shown for volume ratios of 
V c /V = i(+),l(x),2(*),3(H) and 5(H). For all values of the volume ratio V c /V, the pentagonal 
dodecahedron at F = 12 has the same value of S/V 2 / 3 , but for all other F the surface area fluctuates 
widely, although in general it decreases as the volume ratio increases. Also shown is the data for 
bubbles with constant curvature (©) rather than with fixed volume 1 1211 - Shown as horizontal lines 
(from top to bottom) are the value of S/V 2 / 3 for the Kelvin structure (solid line), for the Weaire- 
Phelan structure (quadruple dashes), for the "ideal" flat-faced bubble (triple dashes) and for an 
infinitely large bubble with hexagonal faces (double dashes) lElll . 



4.1 Equal- volume bubble clusters 

We consider first the monodisperse case, relevant to the Kelvin problem, where the volume of the 
central bubble is equal to that of its neighbours, V c = V. As mentioned above, we can go from F = 
5 to 32. The ratio S/V 2 / 3 , shown in figure|3l is lowest at F = 12 (consistent with data for random 
monodisperse foams JHIlX an d increases steeply for F greater than about 16. 

The inset on figure |3] shows the data around the optimal region F = 1 1 to 16. These bubbles, 
which do not pack to fill space, have lower area than Kelvin's (5.306) and even Weaire-Phelan's 



(5.288) (see [22] for details of other space-filling foam structures). They are barely above the value 
for the so-called "ideal" bubble (5.254) [23]. The latter, with F = 13.39, describes a regular (but 
unphysical) "bubble" which would have flat faces, and hence a growth-rate of zero. 

Also of interest is the normalized line-length L/V 1 / 3 of each bubble, plotted in figureHl Note 
that all data lies close to a line L/V 1 / 3 \[F lfl2il . We therefore show the ratio L/V 1 / 3 /F 1 / 2 in 
figure[4| the maximum deviation (i.e. the shallow minimum in the data) occurs for F m 25. 



6 



c 



4.5 



4.45 - 



4.4 - 



4.35 



4.3 



4.25 







— I 1 1 1— 




_L 



X + 
+ 



r 8 







10 20 30 40 50 60 



g5P 



+ ++ ^x x 



X 



X 



xxx^xx ar b. 



*X' 



X 



□ 
□ 



□ 



X* 



9 



^□□■■ B iS° D 



_l I 1_ 



_l I I l_ 



10 



20 30 40 

Number of faces F 



50 



60 



Figure 4: The line-length ratio L/V 1 / 3 increases in proportion to the square-root of F IU2I1 (inset). 
We therefore plot the ratio L/V l / 3 /F 1 / 2 , for volume ratios of V c /V = l(x),2(*),3(H) and 

5(B). The data is everywhere close to L = 4.35V 1/ ' 3 i 7l//2 , confirming the square-root behaviour. 
The dispersity increases as both F and V c /V increase. 



4.2 Non-equal volumes 

4.2.1 Simple volume ratios 

We next consider the case where the volume of the central bubble is not equal to the volume of its 
neighbours. There are again limits to the possible stable clusters, but they vary with the volume 
ratio. We study the simple ratios V c /V = j,2,3 and 5. This choice of volume ratios allows us to 
explore F from 4 to 60. 

Note that the possible range of F is not always continuous. For instance, we cannot construct 
a stable cluster with 26 neighbours for V c /V = j, hence we find F 6 [4 — 25,27]. Similarly, 11 
neighbours is unstable for V c /V = 5, and we find F G [10, 12 — 60]. 

For each value of F we record the topology of each bubble, collated for all volume ratios (table 
[l])> using the notation n x to mean that the bubble has x-faces with n-sides. The topology of the 
central bubble might depend on V c /V: we find such non-uniqueness in only two instances. We 
accept this as due to the slight randomness used in placing the 3F outer points to truncate the 
initial Voronoi pattern. 

The line-length, shown in figure fall close to the same curve as in the monodisperse case. 
The square-root approximation becomes slightly worse as the bubbles become larger and gain 
more faces, with the maximum deviation occurring at higher F for increasing V c /V. 

4.2.2 Large volume ratios 

With larger V c /V we can look at bubbles with many faces and very low surface areas. For instance, 
with F = 122 (corresponding to the fullerene C240) and V c /V = 200 we find a bubble with topology 
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Table 1: The topology of each central bubble, where n x denotes the number* of n- sided faces. 
* denotes configurations for F = 2 and 3 from Hill . 

** denotes alternative configurations for given F with different V c /V : F = 1 1 (42586i for V c /V = 3) 
and 34 (5i 4 6 2 o for V c /V = 2). 
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5i 2 6no and 5/y 2 / 3 = 5.239: see figureEl 

We could extend this process to larger bubbles with more faces. The normalized area should 
eventually approach the value for an infinitely large bubble with hexagonal faces: S/V 2 ^ = 5.229 



1 21]; note that this is not the theoretical lower bound for the normalized area, which corresponds 
to a spherical bubble with F — 1 infinitesimally small neighbours 

4.2.3 Correlations 

Real foams often have a distribution of bubble volumes, and their topology is correlated to the 
geometry: larger bubbles tend to have more neighbours Il25ll . 

Such correlations appear in our results, although we did not specifically include them. Their 
physical origin is clear. In fact, consider a bubble of volume V c , and consider the average of its 
neighbours' volumes, denoted V (mean-field description). Then, for this given V c /V ratio, the 
physically realizable values of F are limited. Within the possible F, the S(F) curves admit an 
optimum: there is a value of F which minimizes the bubble area. These optimal F values do 
increase with V c /V. Moreover, on figure[3]we can read the optimal surface S opt /V 2 ^ as a function 
of F: it is the envelope of all curves plotted, shown in figure |3 It decreases roughly as one over 
the square-root of F as the volume ratio increases. 

In 2D, the expression for L op t/A 1//2 versus n has been used to estimate the energy of a 2D 
foam H2611 . then to determine the correlations between geometry (area A) and topology (number 
of sides n) |EtTi . Here, its 3D counterpart, the S ov i/V 2 ^ versus F relation, is more complicated 



(in particular, unlike in 2D, it depends on the volume ratio) [21]: but it appears to have the same 
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Figure 5: The optimal normalized surface area S opt /V 2 ^ 3 for a range of values of the volume ratio 
V c /V. Data are shown for values of V c /V = 1, 2, 3 and 5 (+) and for the representative calculation 
for F = 122 with V c /V = 200 (x). The limiting value for V c /V -»• °° at S/V 2 / 3 = 5.229 is shown 
as a horizontal line and we also show a power- law fit: S/V 2 / 3 = 5.229 + 0.078/F 423 . 



essential property as in 2D, namely to be a non-increasing function of F; we thus hope to extend 
to 3D this 2D result H- 

In the theory of foam drainage, in which liquid flows along the edges separating the faces 
(Plateau borders), and the coupling of drainage with coarsening, it is useful to know the following 
two dimensionless parameters iflL l28h : V /P and S// 2 , where / is the average length of an edge 
in an F— faced bubble. We can calculate these quantities from our results, and they are shown in 
figure both increase strongly with the number of faces F and are insensitive to the size of the 
neighbouring bubbles. 



4.3 Growth-rate 

As a result of these simulations, we are able to calculate the instantaneous growth-rate of many 
bubbles, with many different numbers of sides, through the formula ©. It is shown in the inset 
to figure [7]- all data lies close to ©, except at (for us unobtainable) small F where the results of 



Cox and Fortes! 1 12 ] are useful 



More instructive is the difference between the analytic formula and our data, shown in figure 
|7J For F > 12, our data are above and below the analytic line: it agrees with the suggestion 



that the analytic formula approximates the average growth-rate 111 111 , and quantifies the dispersion 
around this average (less than 1 % dispersion). Conversely, for F < 12, our data are clustered and 
significantly (up to 10 %) larger than the analytical formula, which confirms that the analytical 
approximations gradually lose their validity at low F, as expected [11]. In a coarsening foam, 
the bubbles with low F are important because it is these bubbles that disappear. So although for 
F > 12 the growth-rate is well approximated by ©, we give in tableElthe growth-rates for bubbles 
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Figure 6: The volume V and surface area S of a bubble with F faces scaled with the average length 
of an edge /. The data for the surface area is affine, S/l 2 = 2.573F -9.801, while for low F the 
volume data is approximately quadratic, V/P « 0.053F 2 . Data is shown for all volume ratios 
V c /V considered here, and the inset shows details of V and S for bubbles with 1 1 to 16 faces. 

Table 2: The growth-rates, averaged over all simulations, for bubbles with few faces, F < 12. They 
differ significantly from the analytic equation © 111 ill , but show very little dispersion. 
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with F < 12, averaged over all simulations. These are almost indistinguishable from the growth 
rates calculated exactly on ideal bubbles [21]: the agreement for F = 4 to 1 1 is better than 1%, and 
it is even lower than 0.1% for regular bubbles (F = 4 and 6). 

5 Conclusions 

The structure of a foam in equilibrium minimizes its (free) energy, which is the product of (i) two 
quantities characterising the system (surface tension and average surface area) and (ii) some func- 
tion of shape only. The structure changes due to coarsening. The coarsening rate is the product of 
a diffusion constant (which depends on the material parameters, including chemical composition), 
that sets the characteristic time scale, and a function only of geometry. Here, we don't consider 
two other phenomena, drainage and film breakage, which cause deviations from equilibrium. 

Using the Surface Evolver, we have studied finite clusters of bubbles to give information about 
the structure of three-dimensional foam and a 3D coarsening law. This approach allows us to get 
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Figure 7: The difference in the rate of change of volume of a bubble with F faces, calculated 
from our simulations using <[T]), and the value calculated according to the analytic formula ©. The 
inset shows the values themselves, again next to the analytic line, from which it deviates at small 
F. Data is shown for volume ratios of V c /V = l(x),2(*),3(H) and 5(B). The data for 

bubbles with constant curvature (©) 1 1211 - rather than fixed volume, is more scattered, but useful 
for low F. 
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a high level of detail and accuracy of the relevant quantities (surface area, pressure difference) 
to get a good insight into how foams coarsen. Our calculated values of the growth law require no 
assumption about the curvature being small, and can be found for bubbles with an arbitrary number 
of faces. 

As the volume ratio between the central bubble and its neighbours changes, we find upper 
and lower bounds on the possible number of faces, because the bubbles deform. This leads us 
to conjecture a value for the kissing problem for foams: no more than 32 bubbles can be stably 
packed around one other of the same volume. 

Although we don't tackle infinite (or, equivalently, periodic) structures, we expect that this data 
will eventually lead to greater insight into the Kelvin problem, since we are starting to understand 
better what happens for bubbles with between 12 and 16 faces. 
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